Lab 4: Nonlinear regression, ARIMA models, ARMAX models, Exponential Smoothing

Sys.setlocale('LC_ALL', 'C')

[1] "LC_CTYPE=C;LC_NUMERIC=C;LC_TIME=C;LC_COLLATE=C;LC_MONETARY=C;LC_MESSAGES=en_US.UTF-8;LC_PAPER=vi_VN;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=vi_VN;LC_IDENTIFICATION=C"

library(fpp2)
library(forecast) 
library(readxl)
library(DIMORA)
setwd('/home/hieu/BDMA/TSA')

Competition between cassettes and CDs

music<- read_excel('music.xlsx')

colors <- grey(seq(0, 1, length = 5))

year<- music$year[1:36]
cass<- music$cassette[1:36]
cd<- music$cd[10:36]
#
##Exploratory plots
#instantaneous
plot(year[10:36],cass[10:36],xlab='Year', ylab='Annual sales',   type= 'b', pch=16, lty=3, cex=0.7, ylim=c(0,1000), col=colors[1])
points(year[10:36],cd, type= 'b', pch=2, lty=3, cex=0.7, col=colors[3])

legend('topleft', legend=c('Cassette', 'CD'), pch=c(16,2), col=c(colors[1], colors[3]))

#
#cumulative
plot(year[10:36],sum(cass[1:9])+cumsum(cass[10:36]),xlab='Year', ylab='Cumulative sales',   type= 'b', pch=16, lty=3, cex=0.7, ylim=c(0,15000), col=colors[1])
points(year[10:36],cumsum(cd), type= 'b', pch=2, lty=3, cex=0.7, col=colors[3])

legend('topleft', legend=c('Cassette', 'CD'), pch=c(16,2), col=c(colors[1], colors[3]))

# UCRCD model with delta and gamma (par = 'double')
ucrcd<- UCRCD(cass,cd, display=T)

summary(ucrcd)

Call: ( UCRCD Model )

  UCRCD(series1 = cass, series2 = cd, display = T)

Residuals Series 1:
    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-61.911  -5.930   2.142   1.660  14.593  55.344 

Residuals Series 2:
     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-84.0608 -13.5927   0.8525   0.1994  18.4845  57.6658 

Coefficients:
 

---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Residual standard error Series 1: 28.17771  on  26  degrees of freedom
 Residual standard error Series 2: 43.43298  on  17  degrees of freedom
 Multiple R-squared:   0.989805  Residual sum of squares:  52712.77

# Plot observed vs fitted
plot(year[10:36],cass[10:36],xlab='Year', ylab='Annual sales',   type= 'b', pch=16, lty=3, cex=0.7, ylim=c(0,1000), col=colors[1])
points(year[10:36],cd, type= 'b', pch=2, lty=3, cex=0.7, col=colors[3])

lines(year[10:36], (ucrcd$fitted.i[[1]])[10:36], lwd=1, col=colors[1])
lines(year[10:36], (ucrcd$fitted.i[[2]]), lwd=1, col=colors[3])

legend('topleft', legend=c('Cassette', 'CD'), lty=1, lwd=2, col=colors[c(1,3)], pch=c(16,2), cex=0.7)

plot(year[10:36],sum(cass[1:9])+cumsum(cass[10:36]),xlab='Year', ylab='Cumulative sales',   type= 'b', pch=16, lty=3, cex=0.7, ylim=c(0,15000), col=colors[1])
points(year[10:36],cumsum(cd), type= 'b', pch=2, lty=3, cex=0.7, col=colors[3])

lines(year[10:36], cumsum(ucrcd$fitted.i[[1]])[10:36], lwd=1, col=colors[1])
lines(year[10:36], cumsum(ucrcd$fitted.i[[2]]), lwd=1, col=colors[3])

legend('topleft', legend=c('Cassette', 'CD'), lty=1, lwd=2, col=colors[c(1,3)], pch=c(16,2), cex=0.7)

ARMAX models

Data on US personal consumption and income

uschange

        Consumption      Income  Production      Savings Unemployment
1970 Q1  0.61598622  0.97226104 -2.45270031   4.81031150          0.9
1970 Q2  0.46037569  1.16908472 -0.55152509   7.28799234          0.5
1970 Q3  0.87679142  1.55327055 -0.35870786   7.28901306          0.5
1970 Q4 -0.27424514 -0.25527238 -2.18545486   0.98522964          0.7
1971 Q1  1.89737076  1.98715363  1.90973412   3.65777061         -0.1
1971 Q2  0.91199291  1.44733417  0.90153584   6.05134180         -0.1
1971 Q3  0.79453885  0.53181193  0.30801942  -0.44583221          0.1
1971 Q4  1.64858747  1.16012514  2.29130441  -1.53087186          0.0
1972 Q1  1.31372218  0.45701150  4.14957387  -4.35859438         -0.2
1972 Q2  1.89147495  1.01662441  1.89062398  -5.05452579         -0.1
1972 Q3  1.53071400  1.90410126  1.27335290   5.80995904         -0.2
1972 Q4  2.31829471  3.89025866  3.43689207  16.04471706         -0.3
1973 Q1  1.81073916  0.70825266  2.79907636  -5.34886849         -0.3
1973 Q2 -0.04173996  0.79430954  0.81768862   8.42603436          0.0
1973 Q3  0.35423556  0.43381827  0.86899693   2.75879565         -0.1
1973 Q4 -0.29163216  1.09380979  1.47296187  11.14642986          0.1
1974 Q1 -0.87702794 -1.66168482 -0.88248358  -2.53351449          0.2
1974 Q2  0.35113555 -0.93835321  0.07427919  -6.59264464          0.3
1974 Q3  0.40959770  0.09448779 -0.41314971   0.51717884          0.5
1974 Q4 -1.47580863 -0.12259599 -4.06411893  11.34339540          1.3
1975 Q1  0.83225762 -0.16369546 -6.85103912  -5.47619069          1.4
1975 Q2  1.65583461  4.53650956 -1.33129558  24.30960536          0.2
1975 Q3  1.41942029 -1.46376532  2.42435972 -17.65616104         -0.4
1975 Q4  1.05437932  0.76166351  2.16904208   0.64809041         -0.2
1976 Q1  1.97998024  1.16825761  3.02720471  -2.95006644         -0.6
1976 Q2  0.91391607  0.51729906  1.27881101  -1.47455755          0.0
1976 Q3  1.05532326  0.73370026  1.30386487  -0.06754475          0.0
1976 Q4  1.29889825  0.59458339  1.77537765  -3.57672239          0.2
1977 Q1  1.13637586 -0.03108003  2.05516067  -9.16055658         -0.4
1977 Q2  0.54994073  1.23808955  3.05838507   9.09050404         -0.2
1977 Q3  0.94985262  1.51880293  1.10308888   7.94495719         -0.4
1977 Q4  1.49599724  1.91456240  0.63346850   6.69627648         -0.4
1978 Q1  0.57549599  0.70266687 -0.29339056   2.92296383         -0.1
1978 Q2  2.11120960  0.98314132  3.94815264  -6.81114259         -0.4
1978 Q3  0.41796279  0.71992620  0.87114701   4.79207162          0.1
1978 Q4  0.79792710  0.78553605  1.78447991   2.37118400          0.0
1979 Q1  0.50584598  1.05755946  0.42594327   7.77418337         -0.2
1979 Q2 -0.05775339 -0.86765105 -0.20491944  -5.28634896         -0.1
1979 Q3  0.97730010  0.47100340 -0.29723637  -1.84549644          0.2
1979 Q4  0.26826982  0.44037974  0.33560928   4.04959810          0.1
1980 Q1 -0.15391875  0.33827686  0.41056141   5.86168864          0.3
1980 Q2 -2.27411019 -1.46388507 -4.30076832   8.24322919          1.3
1980 Q3  1.07188123  1.21301507 -1.64181977   5.70775044         -0.1
1980 Q4  1.31644941  1.94243865  3.78045520   9.15098787         -0.3
1981 Q1  0.52472770 -0.26813406  0.24627687  -5.68139002          0.2
1981 Q2 -0.01728203 -0.02363025  0.30977573   0.88183993          0.1
1981 Q3  0.40165150  2.02680183  0.91707444  15.99035721          0.1
1981 Q4 -0.75287620  0.19560628 -2.25457797   7.80550650          0.9
1982 Q1  0.65938376  0.11969888 -2.07131293  -3.34243955          0.5
1982 Q2  0.36854173  0.57548997 -1.24766384   2.19400166          0.6
1982 Q3  0.76954464  0.53484410 -1.40050430   0.03499563          0.5
1982 Q4  1.80876006  0.44938311 -1.90375664  -9.57651468          0.7
1983 Q1  0.96802954  0.85588425  1.14655720   0.34595460         -0.5
1983 Q2  1.95946831  0.70632719  2.17942248 -10.17004699         -0.2
1983 Q3  1.73949442  1.49810999  3.36771897   0.21217916         -0.9
1983 Q4  1.56389332  2.13138911  2.58168445   8.21600068         -0.9
1984 Q1  0.84526442  2.02348788  2.89709545  13.86918150         -0.5
1984 Q2  1.41504495  1.64921136  1.53821324   4.38900229         -0.6
1984 Q3  0.76546608  1.36163845  0.72128740   6.51686089          0.1
1984 Q4  1.31380062  0.81927319  0.04115557  -2.87544931          0.0
1985 Q1  1.68655320 -0.23895759  0.32353159 -18.71008389         -0.1
1985 Q2  0.93436990  1.90677905  0.07020996  11.82871950          0.2
1985 Q3  1.90256675 -0.33536283 -0.14046924 -23.57393474         -0.3
1985 Q4  0.25656565  1.14181151  0.57978813  11.36628338         -0.1
1986 Q1  0.84304279  1.23951110  0.58132135   5.86126836          0.2
1986 Q2  1.11177390  1.31938549 -0.57641778   3.27551734          0.0
1986 Q3  1.79499406  0.70477150  0.37249329 -10.09044542         -0.2
1986 Q4  0.63768446  0.17977925  1.13734778  -4.82920131         -0.4
1987 Q1  0.01569397  0.81973366  1.30758228  12.46424452          0.0
1987 Q2  1.37731686 -0.97505791  1.75000563 -29.52866718         -0.4
1987 Q3  1.15225712  1.80185055  1.84366200  12.32810406         -0.3
1987 Q4  0.21016439  1.32743427  2.40645058  16.63076101         -0.2
1988 Q1  1.76316026  1.44861875  0.92013121  -0.96896505          0.0
1988 Q2  0.73053714  1.02084894  0.87316353   5.67776867         -0.3
1988 Q3  0.85083233  0.95820336  0.38103668   3.64649867          0.0
1988 Q4  1.13789838  0.96207024  0.70292025  -0.19730358         -0.1
1989 Q1  0.46064152  1.22693023  0.43372685  10.01461545         -0.3
1989 Q2  0.46937808 -0.29489091 -0.36675732  -8.15576525          0.3
1989 Q3  0.98950145  0.67822897 -0.62142121  -2.48622554          0.0
1989 Q4  0.43942767  0.80025832  0.42443392   5.44681102          0.1
1990 Q1  0.85543417  0.83939484  0.68265169   2.87544931         -0.2
1990 Q2  0.31230451  0.59572848  0.77446547   5.10951644          0.0
1990 Q3  0.40261313  0.03740765  0.41944800  -3.17767248          0.7
1990 Q4 -0.75910716 -0.79479735 -1.57345296  -0.17953326          0.4
1991 Q1 -0.34535008  0.21183290 -1.91422028   6.49315257          0.5
1991 Q2  0.83564224  0.69043356  0.59131506  -0.30920615          0.1
1991 Q3  0.48439843  0.36205181  1.36255645  -0.14086493          0.0
1991 Q4 -0.02626579  0.85100324  0.21710308  11.34193010          0.4
1992 Q1  1.85996999  2.12421067 -0.13365365   7.23265150          0.1
1992 Q2  0.68354371  1.04095059  1.76874773   5.46708666          0.4
1992 Q3  1.07661214  0.43562041  0.76167388  -5.93646090         -0.2
1992 Q4  1.18372396  0.34210852  1.05024577  -5.88618856         -0.2
1993 Q1  0.37817936  0.55877186  0.87901471   2.63464703         -0.4
1993 Q2  0.89392729  0.17627103  0.21755108  -6.91664675          0.0
1993 Q3  1.09813766  0.05868803  0.40135891 -11.99337844         -0.3
1993 Q4  0.88122025  0.65496353  1.49618275  -1.83708870         -0.2
1994 Q1  1.14064791  0.69846579  1.22213656  -5.18600629          0.0
1994 Q2  0.77176225  1.05367166  1.78250275   5.15609751         -0.4
1994 Q3  0.77214364  0.59247377  1.26718100  -2.42215898         -0.2
1994 Q4  1.07014805  1.38110661  2.04370404   6.32351898         -0.4
1995 Q1  0.26420505  0.94873528  1.02552601  10.11514398         -0.1
1995 Q2  0.89311141  0.22780635  0.33785685 -10.60541172          0.2
1995 Q3  0.91264702  0.88957006  0.90043887  -0.11570727          0.0
1995 Q4  0.70025425  0.57591998  0.87467273  -2.90726686          0.0
1996 Q1  0.92360967  0.95255663  0.69285195   2.55933958         -0.1
1996 Q2  1.07997887  0.95161791  2.11134752  -0.75802112         -0.2
1996 Q3  0.60055799  0.79369738  1.24418680   3.33843952         -0.1
1996 Q4  0.78298122  0.52035746  1.35396890  -3.33843952          0.2
1997 Q1  1.04949253  0.99858552  1.86714700   0.61269338         -0.2
1997 Q2  0.45219855  0.85103564  1.48763922   6.17532322         -0.2
1997 Q3  1.69654264  1.18352222  2.28632066  -7.22796452         -0.1
1997 Q4  1.18062797  1.42325742  2.48091341   5.43456565         -0.2
1998 Q1  1.02693626  2.10753052  1.10343775  19.35335228          0.0
1998 Q2  1.75069399  1.38767133  0.65122238  -4.81709478         -0.2
1998 Q3  1.30596977  1.01464427  0.72551955  -3.12983982          0.1
1998 Q4  1.45888615  0.80893032  1.44421674  -9.14923404         -0.2
1999 Q1  0.94821191  0.89173174  1.10341663   1.88735718         -0.2
1999 Q2  1.46971415  0.24722185  0.98574261 -23.49652903          0.1
1999 Q3  1.12921436  0.66729226  0.90279881  -9.86264835         -0.1
1999 Q4  1.45748895  1.46092242  1.75533234   2.35825225         -0.2
2000 Q1  1.51106759  1.95061335  0.99682019  12.28684080          0.0
2000 Q2  0.95508878  1.03174349  1.23293805   1.28001748          0.0
2000 Q3  0.96797647  1.16178668 -0.10225268   2.57390229         -0.1
2000 Q4  0.88629738  0.33725343 -0.20388383 -13.16296208          0.0
2001 Q1  0.42159086  0.84865826 -1.35143911  13.22491995          0.4
2001 Q2  0.25689982 -0.08818148 -1.25954437  -6.89043916          0.2
2001 Q3  0.36381084  2.33678920 -1.44101744  41.66826457          0.5
2001 Q4  1.51630321 -1.24443353 -1.06013675 -56.75209674          0.7
2002 Q1  0.29958257  2.40331419  0.70916406  50.75796205          0.0
2002 Q2  0.50899032  0.50559877  1.54280957   0.87861837          0.1
2002 Q3  0.69667241 -0.12828194  0.59478143 -14.70397426         -0.1
2002 Q4  0.53634306  0.47941927 -0.05776556   1.58733492          0.3
2003 Q1  0.43826169  0.27834026  0.53922789   0.49744834         -0.1
2003 Q2  1.10719086  1.43729445 -0.69876172   7.00891625          0.4
2003 Q3  1.46377882  1.62544947  0.60727351   6.18413150         -0.2
2003 Q4  0.77334046  0.40353864  1.00599126  -6.89274778         -0.4
2004 Q1  0.96768535  0.72653162  0.65792806  -2.96152040          0.1
2004 Q2  0.64760607  0.98056746  0.57461780   8.30885627         -0.2
2004 Q3  0.95117167  0.52450113  0.56330030  -8.99318286         -0.2
2004 Q4  1.02041702  1.24238706  1.38522763   6.23585017          0.0
2005 Q1  0.76172556 -0.96827007  1.39435718 -42.28191228         -0.2
2005 Q2  1.08136588  0.78835467  0.50586367 -18.27592893         -0.2
2005 Q3  0.77186494  0.51136949 -0.50305848  -7.87665229          0.0
2005 Q4  0.37591485  0.82191843  0.93365010  20.37236078         -0.1
2006 Q1  1.11522822  2.25904474  0.95057853  37.40653542         -0.2
2006 Q2  0.53100554  0.14987813  0.59636010 -12.34810568         -0.1
2006 Q3  0.58208747  0.28490722  0.33552773 -10.55276140         -0.1
2006 Q4  1.01434389  1.30059162  0.25603401   6.03100080         -0.1
2007 Q1  0.52486184  0.65373993  0.91794957   6.60516929          0.0
2007 Q2  0.33874119  0.19260870  1.19594247  -7.23648452          0.2
2007 Q3  0.44391875  0.26238732  0.22356909  -9.00674555          0.1
2007 Q4  0.12505584  0.08392938  0.16424632   2.32887238          0.3
2008 Q1 -0.20652548  0.71926565 -0.42872571  29.83728599          0.1
2008 Q2  0.16783443  2.08693775 -1.41297022  46.43989041          0.5
2008 Q3 -0.72499446 -2.32611860 -3.26349945 -32.53252494          0.5
2008 Q4 -1.21068558  0.64019534 -4.35417741  36.31240490          1.2
2009 Q1 -0.34354370 -0.18888849 -5.75045075   0.92306020          1.4
2009 Q2 -0.45174364  0.70899368 -3.00372447  16.09059408          0.8
2009 Q3  0.60491332 -1.10343180  1.39880419 -24.49229966          0.3
2009 Q4 -0.01115014 -0.13213193  1.54400617   0.84829220          0.1
2010 Q1  0.53481740  0.10094986  1.88006931  -5.54399051          0.0
2010 Q2  0.81040406  1.29229259  2.05402479  11.65612884         -0.5
2010 Q3  0.64501881  0.49678098  1.42683671  -0.35208609          0.1
2010 Q4  1.01833874  0.69495229  0.37927209  -3.27335958         -0.2
2011 Q1  0.50041315  1.21571502  0.50174040  14.33860193         -0.3
2011 Q2  0.20141978 -0.15658108  0.21878696  -4.07705131          0.1
2011 Q3  0.43372599  0.52891255  1.01113866   2.72250400         -0.1
2011 Q4  0.33593895  0.06074719  0.85151692  -3.45447712         -0.5
2012 Q1  0.60108995  1.62204885  0.88651817  17.62530510         -0.3
2012 Q2  0.16942956  0.76689543  0.62923586   8.96949710          0.0
2012 Q3  0.26416034 -0.05071452  0.07880166  -3.04922177         -0.4
2012 Q4  0.27877186  2.59106697  0.63305509  29.04670355          0.1
2013 Q1  0.46861292 -4.26525047  0.67713243 -68.78826698         -0.4
2013 Q2  0.20545802  0.58146541  0.30744961   7.81647729          0.0
2013 Q3  0.46641787  0.58328912  0.23440888   3.49400682         -0.3
2013 Q4  0.83917367  0.21494896  0.79208722 -11.27661450         -0.5
2014 Q1  0.47345118  1.10369487  0.54709166  13.52020248          0.0
2014 Q2  0.93375698  1.29390492  1.33801074   8.24404770         -0.6
2014 Q3  0.91687178  0.99853396  0.62352731   2.46195256         -0.2
2014 Q4  1.12533250  1.04641801  0.90355427  -1.51305022         -0.3
2015 Q1  0.59624005  0.49040680 -0.46710878  -0.75840017         -0.2
2015 Q2  0.70814389  0.95495949 -0.69702162   5.02391773         -0.1
2015 Q3  0.66496956  0.80166267  0.38060610   3.18092976         -0.3
2015 Q4  0.56167978  0.74006260 -0.84554638   3.48278601          0.0
2016 Q1  0.40468216  0.51902540 -0.41793048   2.23653405          0.0
2016 Q2  1.04770741  0.72372078 -0.20331883  -2.72150106         -0.1
2016 Q3  0.72959779  0.64470081  0.47491844  -0.57285793          0.0

plot(uschange) # 'classical view'

autoplot(uschange) # 'all the series together'

# if we want to see just the first two series
par(mfrow=c(2,1))
plot(uschange[,1])
plot(uschange[,2])

tsdisplay(uschange[,1])

# to go back to 1 panel 
par(mfrow=c(1,1))

##estimate an ARMAX model
armax1<- Arima(uschange[,1], xreg=uschange[,2], order=c(1,0,1))
res1<- residuals(armax1)
Acf(res1)

fitted(armax1)

            Qtr1        Qtr2        Qtr3        Qtr4
1970  0.78006603  0.76998241  0.77997918  0.49470889
1971  0.74691458  0.95367578  0.76138513  0.87286648
1972  0.90802836  1.05689010  1.35176090  1.67290915
1973  1.16367516  1.23640518  0.78531105  0.79005897
1974  0.05003515  0.02887723  0.37542403  0.39653396
1975  0.01085034  1.17153692  0.24080247  0.93242178
1976  0.99640757  1.06214642  0.99472577  0.93136100
1977  0.85572449  1.09847602  0.97310973  1.02294306
1978  0.89498278  0.84279905  1.06793243  0.86730000
1979  0.87840579  0.42871159  0.56972244  0.68221667
1980  0.57340521  0.10300465  0.11068218  0.60940787
1981  0.43861630  0.52865939  0.79754365  0.41344059
1982  0.18783151  0.46954643  0.49131738  0.58362522
1983  0.95498317  0.89199407  1.24400153  1.39600580
1984  1.33190327  1.08208133  1.06911496  0.85908042
1985  0.74756500  1.31211568  0.74540820  1.23350168
1986  0.94718691  0.91606619  0.83474392  0.93429958
1987  0.92479494  0.36025112  1.11010775  0.99798195
1988  0.81430804  0.96157003  0.86476627  0.84641084
1989  0.94794685  0.53818285  0.70483443  0.79746028
1990  0.71661438  0.71355530  0.52369366  0.36479282
1991  0.31875176  0.32696302  0.46657727  0.60408691
1992  0.72976714  0.84119795  0.68756805  0.75856897
1993  0.87592653  0.66138154  0.68859617  0.87554555
1994  0.85504579  0.95994479  0.80156428  0.92224067
1995  0.86455278  0.58247944  0.78777694  0.75592063
1996  0.80440335  0.82728814  0.84656787  0.72326657
1997  0.82077487  0.83972711  0.79895963  1.04971381
1998  1.16710358  0.96916585  1.05446203  1.01986675
1999  1.07924160  0.86920134  1.03626172  1.13920607
2000  1.24577940  1.07983809  1.01776584  0.81608611
2001  0.89572256  0.59217313  0.96084875  0.19014878
2002  1.19322752  0.61764119  0.49594910  0.67128775
2003  0.61055054  0.79469592  0.91574479  0.81404964
2004  0.83743906  0.89386774  0.73441361  0.90861396
2005  0.51585232  0.87493303  0.84602049  0.85710258
2006  0.99228742  0.64174490  0.64139562  0.81848258
2007  0.75090681  0.61238757  0.57143827  0.52902716
2008  0.57538416  0.68538909 -0.17152326  0.32110519
2009 -0.08758067  0.15775900 -0.18758071  0.29165312
2010  0.32840865  0.65512112  0.58502461  0.65935684
2011  0.85147481  0.51911677  0.58684456  0.49275126
2012  0.77004410  0.60517415  0.38933485  0.88819759
2013 -0.45825324  0.69458103  0.59071556  0.52254602
2014  0.78259775  0.75487573  0.76063918  0.81156876
2015  0.77871811  0.80690024  0.75281487  0.72333875
2016  0.65205361  0.64574770  0.74192959            

plot(uschange[,1])
lines(fitted(armax1), col=2)

armax2<- Arima(uschange[,1], xreg=uschange[,2], order=c(1,0,2))
res2<- residuals(armax2)
Acf(res2)

fitted(armax2)

            Qtr1        Qtr2        Qtr3        Qtr4
1970  0.78312833  0.78018625  0.80479451  0.41930867
1971  0.84727727  0.77029485  0.84699950  0.95331366
1972  0.84462865  1.10359869  1.37654140  1.83317035
1973  1.13786583  1.20826089  0.98517520  0.70586675
1974 -0.05888372 -0.10652716  0.15252516  0.37252999
1975  0.26206496  1.00859262  0.13720674  0.91699737
1976  1.21948453  1.08521548  1.14267609  0.94879507
1977  0.77480045  1.08771984  1.08065070  0.98559686
1978  0.77381614  0.89892787  0.91640998  1.05941956
1979  0.89261857  0.38091069  0.53759910  0.54387460
1980  0.62276776  0.12597681  0.28963813  0.22350695
1981  0.29439649  0.66462244  1.02513172  0.44125915
1982  0.22393882  0.25312328  0.48707851  0.59336917
1983  0.90313806  1.08127964  1.25264967  1.50432058
1984  1.44036881  1.17436646  0.95189308  0.84432657
1985  0.62302160  1.25286288  0.88979506  1.13370425
1986  1.15901792  0.86946926  0.70947373  0.78744003
1987  1.06836399  0.45881327  0.90243930  1.03659542
1988  0.96270927  0.78370422  0.93047644  0.86924854
1989  0.91533693  0.58627680  0.66536135  0.72658394
1990  0.77608713  0.67685037  0.55014995  0.30925675
1991  0.39935061  0.27248164  0.26732620  0.62816783
1992  0.89742511  0.70072700  0.80056755  0.77549795
1993  0.88877001  0.77843395  0.62254092  0.81818587
1994  0.90212475  0.96475852  0.84927848  0.92193247
1995  0.82302821  0.63143690  0.68951971  0.73146966
1996  0.84367658  0.83076836  0.83169021  0.77671808
1997  0.80639086  0.80452612  0.86904546  0.93489503
1998  1.25190643  1.05908067  0.96434947  1.05049178
1999  1.09414605  0.93073127  0.96307717  1.17990917
2000  1.25691761  1.07448753  1.05545026  0.78640224
2001  0.84779060  0.60439534  0.96407640  0.13244226
2002  0.98024037  0.82330683  0.48567225  0.59802581
2003  0.62215476  0.83382918  0.88479540  0.77211433
2004  0.92453660  0.92649952  0.77051709  0.86179287
2005  0.46441540  0.86776267  0.86193275  0.90781283
2006  1.07577457  0.54136888  0.60248018  0.82077639
2007  0.72101375  0.64631554  0.58322422  0.49104710
2008  0.58814023  0.74866969 -0.28741575  0.27567649
2009 -0.01283404  0.01246222 -0.24998664  0.14688615
2010  0.47625237  0.73744143  0.63427916  0.71035984
2011  0.86308493  0.59824770  0.60961002  0.44522796
2012  0.76540887  0.60120651  0.39972191  0.89227477
2013 -0.52514212  0.53160833  0.74818031  0.56648441
2014  0.74603201  0.83062237  0.73814649  0.80864084
2015  0.76900278  0.88423189  0.77020729  0.70789375
2016  0.63785994  0.64367730  0.65987214            

plot(uschange[,1])
lines(fitted(armax2), col=2)

AIC(armax2)

[1] 325.908

AIC(arima1)

[1] 353.1935

# Note: forecasts are available only if we have future values of 'income'

# procedure also available with auto.arima
auto.arima<- auto.arima(uschange[,1], xreg=uschange[,2])

Quarterly international arrivals (in thousands) to Australia from Japan, New Zealand, the UK and the US. 1981Q1 - 2012Q3

?arrivals
autoplot(arrivals)

autoplot(arrivals[,c(1,2)])

Japan<- arrivals[,1]
NZ<- arrivals[,2]
UK<- arrivals[,3]
US<- arrivals[,4]

# #we try with a simple arima model, ARMAX
auto.a<- auto.arima(NZ, xreg=Japan) 
AIC(auto.a)

[1] 973.6105

# We try with a regression model with trend, season, and external variable 'Japan'
mod<- tslm(NZ~ trend+season+Japan) 
summary(mod)

Call:
tslm(formula = NZ ~ trend + season + Japan)

Residuals:
   Min     1Q Median     3Q    Max 
-40.32 -14.12  -2.51  12.70  55.63 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.73757    5.21347   0.909  0.36531    
trend        2.21843    0.05539  40.052  < 2e-16 ***
season2     36.35773    5.01769   7.246 4.38e-11 ***
season3     63.33026    4.92212  12.866  < 2e-16 ***
season4     45.09142    4.95892   9.093 2.33e-15 ***
Japan       -0.10039    0.03232  -3.106  0.00236 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 19.67 on 121 degrees of freedom
Multiple R-squared:  0.9482,    Adjusted R-squared:  0.946 
F-statistic: 442.9 on 5 and 121 DF,  p-value: < 2.2e-16

fitted(mod)

           Qtr1       Qtr2       Qtr3       Qtr4
1981   5.473922  44.596412  73.702533  56.744161
1982  14.111305  53.340009  82.418522  65.472097
1983  22.628218  62.051280  91.008119  74.004269
1984  31.251144  70.517295  99.486582  82.318292
1985  39.414981  78.991039 108.118945  90.566457
1986  47.510752  86.651414 116.267522  98.337764
1987  55.004376  94.934013 122.563578 104.736118
1988  61.010502  99.657338 127.933523 110.404225
1989  68.491577 108.845668 136.838850 120.604200
1990  74.807108 114.657135 141.688969 126.039901
1991  83.682418 122.719275 149.358178 132.036088
1992  87.828292 128.898776 156.225560 140.169807
1993  95.524004 136.766759 164.328660 147.883890
1994 103.102559 145.007697 171.680632 155.159866
1995 111.462962 152.265100 177.464491 163.072122
1996 117.239090 161.136997 185.808128 172.519362
1997 126.483740 170.719062 193.849186 181.070508
1998 137.181756 179.467177 205.411471 191.863193
1999 147.224618 189.712227 215.568177 201.293668
2000 156.278926 198.601893 225.039210 208.014379
2001 164.989795 207.314770 233.437962 222.448459
2002 175.508110 215.982270 242.616856 225.373172
2003 184.495357 230.450182 254.243793 234.589712
2004 192.461222 234.231033 260.581310 243.709174
2005 200.494450 244.901943 269.998834 253.592614
2006 209.805860 254.432307 280.136766 263.547634
2007 221.025108 265.056434 290.474578 274.694299
2008 232.838170 276.199485 302.070293 287.263910
2009 244.765880 287.395341 314.315139 297.607746
2010 253.343530 295.716391 320.226796 306.002282
2011 263.841867 306.382683 331.249779 315.888633
2012 272.184400 314.617597 339.578056           

plot(NZ)
lines(fitted(mod), col=2)

plot(residuals(mod))

# analysis of residuals: autocorrelation? 
tsdisplay(residuals(mod))

# fit an arima model to residuals
aar<- auto.arima(residuals(mod))
fitted(aar)

            Qtr1        Qtr2        Qtr3        Qtr4
1981  43.6224119  42.8277172  12.1263286   5.1327012
1982  39.2568501  30.9438192  -7.0955060  -9.0387482
1983  19.5134000   9.3981650 -17.6940294 -17.4483315
1984  15.6074063  -2.5341788 -18.3538201 -19.9414017
1985   1.4237203  -4.5966746 -32.7036950 -33.3063718
1986   7.9409396  -6.3701574 -23.8781040  -8.2219509
1987  18.2048570   8.7177606 -10.6426164   6.1043757
1988  31.7536169  19.4098627   9.5556372  40.0743379
1989  42.7494458  21.4951003   5.9384598 -11.5850634
1990   9.7508942  -3.2840699 -16.1079337  -4.5609618
1991  -2.7173333   2.6010011 -10.8479635  -5.0272017
1992   8.0078896 -17.3846533 -18.6241736 -19.5752213
1993  -6.4796214  -9.0315699 -16.7577210  -4.0981064
1994  -5.3500015 -20.1121535 -24.2384935 -22.4297611
1995 -21.0249533 -27.2382296 -21.7057467 -12.6737140
1996  -5.7705269  -3.2515084   7.0836384   6.5702005
1997  13.2841639   1.9723677   0.9505590  -5.1161518
1998  11.7656134 -19.7140471  10.7822451  -2.7909927
1999  -1.5343071  -0.9442189  -1.6686489 -13.3979079
2000  -1.5534222 -12.9213959  15.9815234  -3.6481466
2001  10.8415074  12.1928739   4.5535895  -0.2197456
2002 -14.2595975 -13.9064412 -15.7578096 -27.9521114
2003 -11.7220191 -21.1521442 -22.1043363 -11.5919890
2004 -11.4794723  -5.6563040  15.6924561  37.1415765
2005  25.0658788  27.7020028  40.5272783  39.4642352
2006  13.4666275   8.5647731  33.0138668  15.9476795
2007   4.0567564  14.3061907  42.3182107  27.2034387
2008   0.0918427  15.4760632  17.3447321  -0.2983433
2009 -11.0766984  -0.1767021  -0.5943938   5.4166843
2010 -24.2731970  -5.7687911  -1.5124909   6.3192243
2011 -14.2285524  -7.6044283   4.0132082  -2.5252514
2012 -26.7724317  -9.4660669  -9.6347055            

# complete the analysis by summing predictions made with linear model and arma on residuals
plot(NZ)
lines(fitted(mod)+fitted(aar), col=3)

# Notice the difference between the two methods: ARMAX and linear regression+ Arima on residuals

# #Another way of performing the same linear regression
tt<- (1:length(NZ))
seas <- factor(c(rep(1:4,length(NZ)/4),1:3)) #1:3 because there are three observations 'out' 
mod2 <- lm(NZ~ tt+seas+Japan)
summary(mod2)

Call:
lm(formula = NZ ~ tt + seas + Japan)

Residuals:
   Min     1Q Median     3Q    Max 
-40.32 -14.12  -2.51  12.70  55.63 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.73757    5.21347   0.909  0.36531    
tt           2.21843    0.05539  40.052  < 2e-16 ***
seas2       36.35773    5.01769   7.246 4.38e-11 ***
seas3       63.33026    4.92212  12.866  < 2e-16 ***
seas4       45.09142    4.95892   9.093 2.33e-15 ***
Japan       -0.10039    0.03232  -3.106  0.00236 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 19.67 on 121 degrees of freedom
Multiple R-squared:  0.9482,    Adjusted R-squared:  0.946 
F-statistic: 442.9 on 5 and 121 DF,  p-value: < 2.2e-16

AIC(mod2)

[1] 1124.947

AIC(mod)

[1] 1124.947

Forecasting electricity demand

elecdaily

Time Series:
Start = c(1, 4) 
End = c(53, 4) 
Frequency = 7 
            Demand WorkDay Temperature
 1.428571 174.8963       0        26.0
 1.571429 188.5909       1        23.0
 1.714286 188.9169       1        22.2
 1.857143 173.8142       0        20.3
 2.000000 169.5152       0        26.1
 2.142857 195.7288       1        19.6
 2.285714 199.9029       1        20.0
 2.428571 205.3375       1        27.4
 2.571429 228.0782       1        32.4
 2.714286 258.5984       1        34.0
 2.857143 201.7970       0        22.4
 3.000000 187.6298       0        22.5
 3.142857 254.6636       1        30.0
 3.285714 322.2323       1        42.4
 3.428571 343.9934       1        41.5
 3.571429 347.6376       1        43.2
 3.714286 332.9455       1        43.1
 3.857143 219.7517       0        23.7
 4.000000 186.9816       0        22.3
 4.142857 228.4876       1        24.0
 4.285714 222.3650       1        23.1
 4.428571 220.4669       1        23.3
 4.571429 240.5806       1        29.6
 4.714286 224.6652       1        22.1
 4.857143 181.1498       0        20.3
 5.000000 182.2231       0        27.0
 5.142857 230.2226       0        34.5
 5.285714 309.3744       1        41.4
 5.428571 243.2239       1        24.5
 5.571429 258.3896       1        29.2
 5.714286 259.3053       1        26.8
 5.857143 241.3185       0        29.2
 6.000000 264.5050       0        38.6
 6.142857 264.9556       1        33.9
 6.285714 221.6807       1        25.7
 6.428571 237.1498       1        30.5
 6.571429 276.4134       1        35.1
 6.714286 283.3155       1        35.6
 6.857143 272.5084       0        40.0
 7.000000 243.6298       0        39.5
 7.142857 223.0845       1        24.7
 7.285714 238.7871       1        26.2
 7.428571 252.2366       1        26.9
 7.571429 250.9030       1        26.6
 7.714286 242.4558       1        26.6
 7.857143 214.2726       0        28.7
 8.000000 175.2935       0        23.0
 8.142857 216.4104       1        25.3
 8.285714 232.6376       1        30.2
 8.428571 222.3622       1        25.4
 8.571429 210.2031       1        20.3
 8.714286 211.6370       1        20.4
 8.857143 186.7995       0        19.6
 9.000000 178.2322       0        22.9
 9.142857 221.6854       1        24.3
 9.285714 237.5595       1        32.7
 9.428571 220.5093       1        21.8
 9.571429 216.7880       1        21.5
 9.714286 214.8875       1        23.5
 9.857143 193.5388       0        23.1
10.000000 183.3550       0        22.3
10.142857 228.3825       1        25.0
10.285714 255.3741       1        33.3
10.428571 233.9622       1        25.6
10.571429 216.4703       1        20.8
10.714286 215.9607       1        22.5
10.857143 194.5699       0        27.8
11.000000 200.7804       0        32.7
11.142857 213.9160       0        32.1
11.285714 240.5818       1        31.1
11.428571 220.3766       1        21.8
11.571429 219.1226       1        22.7
11.714286 223.9310       1        26.2
11.857143 191.9780       0        29.1
12.000000 172.9602       0        20.3
12.142857 212.6693       1        24.7
12.285714 216.7992       1        24.1
12.428571 215.9562       1        22.2
12.571429 226.8287       1        30.7
12.714286 221.9756       1        22.9
12.857143 186.2491       0        19.4
13.000000 179.2746       0        18.7
13.142857 213.5667       1        18.8
13.285714 214.6133       1        24.4
13.428571 217.4129       1        26.9
13.571429 221.7905       1        25.3
13.714286 214.8814       1        21.5
13.857143 185.2121       0        23.6
14.000000 181.7951       0        25.9
14.142857 230.6637       1        30.1
14.285714 253.2901       1        32.9
14.428571 235.5289       1        26.0
14.571429 222.7004       1        18.2
14.714286 218.9320       1        19.7
14.857143 192.1518       0        24.1
15.000000 182.5951       0        24.3
15.142857 218.4745       1        24.5
15.285714 226.2018       1        22.0
15.428571 225.5993       1        17.8
15.571429 223.8808       1        18.9
15.714286 220.3987       1        18.9
15.857143 189.7302       0        22.4
16.000000 180.7974       0        21.5
16.142857 212.9014       1        20.7
16.285714 216.1034       1        21.1
16.428571 218.3603       1        21.6
16.571429 215.1137       1        24.3
16.714286 177.7474       0        18.8
16.857143 183.7940       0        17.3
17.000000 178.9498       0        18.3
17.142857 181.3320       0        18.1
17.285714 219.8456       1        17.6
17.428571 218.9037       1        19.2
17.571429 217.9708       1        18.3
17.714286 189.6816       0        20.5
17.857143 189.3522       0        19.2
18.000000 186.1761       0        17.2
18.142857 222.7253       1        21.3
18.285714 226.8940       1        22.3
18.428571 227.8383       1        17.6
18.571429 234.1514       1        15.6
18.714286 234.8755       1        13.9
18.857143 205.5938       0        17.0
19.000000 203.6670       0        13.7
19.142857 237.8857       1        16.4
19.285714 240.6153       1        16.0
19.428571 240.3634       1        16.2
19.571429 241.3502       1        18.4
19.714286 236.7497       1        18.3
19.857143 206.3455       0        16.3
20.000000 192.7160       0        19.7
20.142857 229.3544       1        19.9
20.285714 230.0164       1        20.1
20.428571 224.4637       1        20.2
20.571429 216.9592       1        22.4
20.714286 212.6305       1        22.8
20.857143 188.9991       0        20.0
21.000000 181.5075       0        21.4
21.142857 217.0623       1        22.8
21.285714 224.3771       1        21.6
21.428571 221.9410       1        23.1
21.571429 222.4899       1        20.9
21.714286 221.5659       1        20.3
21.857143 195.6669       0        17.1
22.000000 190.3111       0        20.2
22.142857 226.0510       1        20.8
22.285714 226.9175       1        21.7
22.428571 230.2927       1        18.0
22.571429 234.6604       1        15.8
22.714286 231.8486       1        17.7
22.857143 201.0383       0        18.4
23.000000 199.7129       0        14.9
23.142857 232.5876       1        16.4
23.285714 235.0114       1        16.4
23.428571 230.9970       1        18.1
23.571429 232.0875       1        16.7
23.714286 234.6715       1        15.4
23.857143 202.4323       0        16.8
24.000000 196.7197       0        14.9
24.142857 206.5046       0        15.4
24.285714 235.7782       1        18.3
24.428571 241.0994       1        16.2
24.571429 243.6167       1        16.2
24.714286 237.2126       1        16.5
24.857143 212.1366       0        15.2
25.000000 203.3582       0        16.0
25.142857 241.3206       1        17.0
25.285714 244.7568       1        15.8
25.428571 244.5883       1        14.2
25.571429 252.8974       1        12.2
25.714286 239.6948       1        15.1
25.857143 211.0494       0        15.0
26.000000 205.6888       0        15.4
26.142857 245.3576       1        14.4
26.285714 261.7745       1        14.3
26.428571 248.2476       1        15.7
26.571429 243.3896       1        15.6
26.714286 241.5915       1        14.5
26.857143 218.8741       0        12.7
27.000000 220.2953       0        11.0
27.142857 255.0056       1        13.1
27.285714 254.8101       1        13.1
27.428571 245.0865       1        16.6
27.571429 247.8669       1        12.4
27.714286 249.0520       1        12.4
27.857143 218.4067       0        15.2
28.000000 212.9148       0        13.9
28.142857 241.2600       1        15.7
28.285714 242.9725       1        14.3
28.428571 257.5991       1        13.3
28.571429 257.9806       1        14.7
28.714286 247.6547       1        14.4
28.857143 227.0290       0        11.1
29.000000 220.9397       0        14.9
29.142857 258.8287       1        13.1
29.285714 264.1827       1        12.9
29.428571 253.5232       1        14.6
29.571429 259.8050       1        12.1
29.714286 256.2355       1        13.8
29.857143 223.4825       0        13.2
30.000000 215.7748       0        13.2
30.142857 258.8278       1        12.9
30.285714 268.9846       1        12.2
30.428571 267.2355       1        13.4
30.571429 269.3544       1        14.9
30.714286 253.2748       1        15.9
30.857143 223.1568       0        14.9
31.000000 211.2160       0        16.8
31.142857 242.0725       1        16.7
31.285714 246.4046       1        16.5
31.428571 237.8841       1        18.9
31.571429 239.6187       1        19.3
31.714286 263.3148       1         9.8
31.857143 227.5769       0        13.4
32.000000 221.6350       0        14.2
32.142857 256.6197       1        15.2
32.285714 250.4816       1        16.5
32.428571 246.8155       1        15.9
32.571429 252.0098       1        13.6
32.714286 243.5690       1        15.9
32.857143 217.1874       0        14.6
33.000000 212.9136       0        13.6
33.142857 261.0586       1        12.7
33.285714 262.4751       1        12.7
33.428571 253.5156       1        14.5
33.571429 249.8530       1        14.8
33.714286 244.9549       1        16.9
33.857143 221.1463       0        13.3
34.000000 210.5977       0        15.7
34.142857 250.2505       1        13.6
34.285714 248.2650       1        13.3
34.428571 247.3034       1        13.3
34.571429 242.7937       1        17.4
34.714286 231.1652       1        20.6
34.857143 202.5800       0        19.9
35.000000 200.8607       0        16.3
35.142857 239.1268       1        16.1
35.285714 236.9323       1        18.7
35.428571 233.0740       1        19.5
35.571429 233.1221       1        20.0
35.714286 233.8208       1        15.0
35.857143 193.0105       0        21.0
36.000000 189.3291       0        21.0
36.142857 236.2696       1        17.8
36.285714 238.2240       1        16.1
36.428571 241.4517       1        14.8
36.571429 241.7054       1        15.9
36.714286 233.5166       1        16.7
36.857143 198.1050       0        18.9
37.000000 186.0104       0        20.7
37.142857 218.2276       1        20.8
37.285714 222.7274       1        20.7
37.428571 222.0213       1        19.1
37.571429 226.7373       1        17.1
37.714286 231.3683       1        15.9
37.857143 196.5834       0        18.5
38.000000 185.0438       0        20.1
38.142857 225.4641       1        18.0
38.285714 235.6982       1        14.8
38.428571 236.1329       1        16.2
38.571429 241.1509       1        13.7
38.714286 230.3640       1        15.5
38.857143 199.8989       0        16.8
39.000000 187.0570       0        17.2
39.142857 217.4684       1        23.0
39.285714 215.1564       1        25.3
39.428571 215.4661       1        23.5
39.571429 217.7978       1        17.0
39.714286 214.7996       1        19.3
39.857143 184.7269       0        23.3
40.000000 177.3159       0        27.9
40.142857 207.7305       1        22.9
40.285714 218.1746       1        23.8
40.428571 224.6942       1        15.9
40.571429 224.9024       1        17.4
40.714286 212.3656       1        18.1
40.857143 183.7487       0        26.2
41.000000 173.5839       0        19.2
41.142857 221.1614       1        25.0
41.285714 216.7209       1        17.7
41.428571 221.2251       1        18.7
41.571429 217.7267       1        23.1
41.714286 213.9397       1        24.9
41.857143 186.2191       0        21.2
42.000000 181.4456       0        27.5
42.142857 222.1222       1        15.5
42.285714 228.7036       1        16.9
42.428571 225.6833       1        17.1
42.571429 224.5799       1        16.7
42.714286 214.5551       1        18.1
42.857143 191.7391       0        24.2
43.000000 183.9389       0        28.8
43.142857 213.2147       1        20.6
43.285714 222.9096       1        28.2
43.428571 233.0016       1        31.1
43.571429 224.6348       1        22.7
43.714286 227.9454       1        30.4
43.857143 193.2198       0        19.8
44.000000 182.9607       0        24.5
44.142857 212.9383       1        22.5
44.285714 220.9596       1        19.7
44.428571 218.7275       1        22.5
44.571429 218.7459       1        21.8
44.714286 226.7615       1        32.7
44.857143 187.8364       0        20.3
45.000000 186.3019       0        16.6
45.142857 200.5634       1        23.3
45.285714 186.9384       0        28.9
45.428571 206.6377       1        17.8
45.571429 211.9857       1        23.5
45.714286 227.2132       1        32.8
45.857143 205.4283       0        34.6
46.000000 176.0597       0        19.4
46.142857 211.9239       1        20.2
46.285714 215.3201       1        18.8
46.428571 222.2635       1        24.3
46.571429 239.6066       1        30.6
46.714286 223.6820       1        27.7
46.857143 189.5637       0        18.0
47.000000 182.0131       0        18.1
47.142857 208.3546       1        20.1
47.285714 211.8105       1        20.9
47.428571 217.4702       1        29.6
47.571429 228.7774       1        28.1
47.714286 215.2890       1        23.0
47.857143 196.0592       0        30.8
48.000000 187.2593       0        25.0
48.142857 222.5512       1        19.6
48.285714 213.9397       1        21.8
48.428571 218.0656       1        18.7
48.571429 214.9400       1        19.4
48.714286 212.5017       1        20.5
48.857143 193.7848       0        24.3
 [ reached getOption("max.print") -- omitted 32 rows ]

plot(elecdaily)

# fit a quadratic regression model with ARMA errors using the auto.arima function
xreg <- cbind(MaxTemp = elecdaily[, 'Temperature'],
              MaxTempSq = elecdaily[, 'Temperature']^2,
              Workday = elecdaily[, 'WorkDay'])
fit <- auto.arima(elecdaily[, 'Demand'], xreg = xreg)
checkresiduals(fit)

    Ljung-Box test

data:  Residuals from Regression with ARIMA(2,1,2)(2,0,0)[7] errors
Q* = 28.229, df = 8, p-value = 0.0004326

Model df: 6.   Total lags used: 14

# The model has some significant autocorrelation in the residuals

# using the estimated model we forecast 14 days ahead, setting for the next 14 days a temperature at a constant level of 26.
fcast <- forecast(fit,
                  xreg = cbind(MaxTemp=rep(26,14), MaxTempSq=rep(26^2,14),
                               Workday=c(0,1,0,0,1,1,1,1,1,0,0,1,1,1)))

plot(fcast)

Lagged Predictors

plot(insurance, main='Insurance advertising and quotations', xlab='year', xaxt='n')
axis(1, at=c(2002,2003,2004,2005),labels=c('2002','2003','2004','2005'), line=-1, lwd=0, lwd.ticks=0.5)

# we consider 1 lagged predictor
advert<- cbind(insurance[,2], c(NA,insurance[1:39,2]))
colnames(advert)<- paste('AdLag',0:1, sep='')
fit<- auto.arima(insurance[,1], xreg=advert[,1:2],d=0)
summary(fit)

Series: insurance[, 1] 
Regression with ARIMA(3,0,0) errors 

Coefficients:
         ar1      ar2     ar3  intercept  AdLag0  AdLag1
      1.4117  -0.9317  0.3591     2.0393  1.2564  0.1625
s.e.  0.1698   0.2545  0.1592     0.9931  0.0667  0.0591

sigma^2 = 0.2165:  log likelihood = -23.89
AIC=61.78   AICc=65.4   BIC=73.43

Training set error measures:
                     ME      RMSE       MAE        MPE     MAPE      MASE        ACF1
Training set -0.0126932 0.4344493 0.3501642 -0.1590732 2.611605 0.1020281 0.001076967

# we provide the forecast by assuming that advertising is 8 units in each future month
f<-forecast(fit,h=20, xreg=cbind(AdLag0=rep(8,20),AdLag1=c(advert[40,1],rep(8,19))))
plot(f)

SARMAX refinement for diffusion models

Twitter Revenues

twitter<- read_excel('twitter.xlsx')
length(twitter$twitter)

[1] 46

tw<- (twitter$twitter)

# GGM 
GGM_tw<- GGM(tw, prelimestimates=c(4.463368e+04, 0.001, 0.01, 1.923560e-03, 9.142022e-02))

Warning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs produced

summary(GGM_tw)

Call: ( Guseo Guidolin Model )

  GGM(series = tw, prelimestimates = c(44633.68, 0.001, 0.01, 0.00192356, 
    0.09142022))

Residuals:
    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-257.21  -64.74   20.95    9.17   70.88  265.11 

Coefficients:
 

---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Residual standard error  103.2668  on  41  degrees of freedom
 Multiple R-squared:   0.999862  Residual sum of squares:  437225.1

pred_GGM_tw<- predict(GGM_tw, newx=c(1:60))
pred_GGM_tw.inst<- make.instantaneous(pred_GGM_tw)

plot(tw, type= 'b',xlab='Quarter', ylab='Quarterly revenues',  pch=16, lty=3, cex=0.6, xlim=c(1,60))
lines(pred_GGM_tw.inst, lwd=2, col=2)

lines(pred_GGM_tw.inst, lwd=2, col=3)

# Analysis of residuals
res_GGMtw<- residuals(GGM_tw)
acf<- acf(residuals(GGM_tw))

fit_GGMtw<- fitted(GGM_tw) ##they are cumulative values
fit_GGMtw_inst<- make.instantaneous(fit_GGMtw)

SARMAX refinement

# SARMAX model with external covariate 'fit_GGM' 
s2 <- Arima(cumsum(tw), order = c(3,0,1), seasonal=list(order=c(3,0,1), period=4),xreg = fit_GGMtw)
summary(s2)

Series: cumsum(tw) 
Regression with ARIMA(3,0,1)(3,0,1)[4] errors 

Coefficients:
         ar1      ar2     ar3     ma1    sar1    sar2     sar3    sma1  intercept    xreg
      1.8294  -1.2030  0.1683  -1.000  0.5144  0.6488  -0.2487  0.0400    21.0796  0.9983
s.e.  0.1519   0.2605  0.1508   0.027  0.5208  0.2650   0.3104  0.5359    10.4055  0.0010

sigma^2 = 2329:  log likelihood = -242.57
AIC=507.14   AICc=514.91   BIC=527.26

Training set error measures:
                    ME     RMSE      MAE       MPE     MAPE       MASE       ACF1
Training set 0.2324242 42.69616 30.24614 0.6453786 2.042453 0.04979198 0.03877914

pres2 <- make.instantaneous(fitted(s2))

plot(twitter$twitter, type= 'b',xlab='Quarter', ylab='Quarterly revenues',  pch=16, lty=3, xaxt='n', cex=0.6)
axis(1, at=c(1,10,19,28,37,46), labels=twitter$quarter[c(1,10,19,28,37,46)])

lines(fit_GGMtw_inst, lwd=1, lty=2)
lines(pres2, lty=1,lwd=1, col=2)

Interest in Zoom according to searches in Google ####

interest<- read.csv('zoomgoogle.csv')
int<- interest$zoom[90:180]
week <- as.Date(interest$week[90:180])
# we set a 'timing' more refined for graphical purposes only
tfine150 <- seq(1,150,0.01) 

# Exploratory plots
plot(int, type= 'b',xlab='Week', ylab='Weekly Google searches',  pch=16, lty=3, cex=0.6, xaxt='n')
axis(1, at=c(1,18,36,54,72,90), labels=format(week[c(1,18,36,54,72,90)], '%d/%m/%y'))

plot(cumsum(int), type= 'b',xlab='Week', ylab='Cumulative Google searches',  pch=16, lty=3, cex=0.6, xaxt='n')
axis(1, at=c(1,18,36,54,72,90), labels=format(week[c(1,18,36,54,72,90)], '%d/%m/%y'))

# GBM with one exponential shock
gbme1.go<- GBM(int,shock = 'exp',nshock = 1,prelimestimates = c(4.046997e+03, 0.001, 0.1, 30,-01,5))

# Prediction with GBMe1
pred.gbme1go<- predict(gbme1.go, newx=tfine150)
pred.gbme1goi<- make.instantaneous(pred.gbme1go)[-1]

# Plots observed vs predicted
plot(int, type= 'b',xlab='Week', ylab='Weekly Google searches',  pch=16, lty=3, cex=0.6, xaxt='n', col=colors[3])
axis(1, at=c(1,18,36,54,72,90), labels=format(week[c(1,18,36,54,72,90)], '%d/%m/%y'))

lines(tfine150[-1], pred.gbme1goi*100, lwd=2, col=1)

plot(cumsum(int), type= 'b',xlab='Week', ylab='Cumulative Google searches',  pch=16, lty=3, cex=0.6, xaxt='n',col=colors[3])
axis(1, at=c(1,18,36,54,72,90), labels=format(week[c(1,18,36,54,72,90)], '%d/%m/%y'))

lines(tfine150, pred.gbme1go, lwd=2, col=1)

SARMAX refinement

fit.google<- fitted(gbme1.go)
s2 <- Arima(cumsum(int), order = c(1,0,1), seasonal=list(order=c(0,0,1), period=52), xreg = fit.google)
summary(s2)

Series: cumsum(int) 
Regression with ARIMA(1,0,1)(0,0,1)[52] errors 

Coefficients:
         ar1     ma1     sma1  intercept    xreg
      0.9155  0.5564  -0.1156    -0.1915  1.0019
s.e.  0.0402  0.0782   0.2036    10.9422  0.0071

sigma^2 = 16.99:  log likelihood = -257.29
AIC=526.58   AICc=527.58   BIC=541.65

Training set error measures:
                      ME     RMSE     MAE       MPE     MAPE       MASE      ACF1
Training set -0.02828682 4.007078 2.47122 0.1528441 1.082046 0.08839818 0.1306401

pres2 <- make.instantaneous(fitted(s2))

# Plots observed vs predicted with SARMAX refinement
plot(int, type= 'b',xlab='week', ylab='Weekly Google searches',  pch=16, lty=3, cex=0.6, xaxt='n', col=colors[3])
axis(1, at=c(1,18,36,54,72,90), labels=format(week[c(1,18,36,54,72,90)], '%d/%m/%y'))

lines(tfine150[-1], pred.gbme1goi*100, lwd=1, lty=2)
lines(pres2, lty=1,lwd=1)

legend('topright', legend=c('GBMe1','GBMe1+SARMAX'), lty=c(2,1))

plot(cumsum(int), type= 'b',xlab='week', ylab='Cumulative Google searches',  pch=16, lty=3, cex=0.6, xaxt='n', col=colors[3])
axis(1, at=c(1,18,36,54,72,90), labels=format(week[c(1,18,36,54,72,90)], '%d/%m/%y'))

lines(tfine150, pred.gbme1go, lwd=1, lty=2)
lines(cumsum(pres2), lty=1,lwd=1)

legend('bottomright', legend=c('GBMe1','GBMe1+SARMAX'), lty=c(2,1))